!
!     PHOTOMETRY   This subroutine derives the aperture photometry.
!     Copyright (C) 1997-9  Filip Hroch, Masaryk University, Brno, CZ
!                    
!
!===========================================================================
!
subroutine  photsb (ncol, nrow, d, watch, coofil, magfil)
!
!=======================================================================
!
! This subroutine derives the concentric aperture photometry.  At
! present, this is the only place in all of DAOPHOT where sky values
! are derived for the individual stars.
!
!               OFFICIAL DAO VERSION:  1991 April 18
!
! Argument
!
!    WATCH (INPUT) governs whether information relating to the progress
!          of the reductions is to be typed on the terminal screen
!          during execution.
!
!=======================================================================
!
  implicit none
  integer :: ncol, nrow
  real :: d(ncol, nrow)
  integer, parameter  :: minsky = 20, maxap = 12
  real, parameter :: Pi = 3.14159
  character(len=*) :: coofil, magfil
!
! Parameters:
!
! MINSKY is the smallest number of pixels from which the sky may be
!        determined.  If for some star the number of sky pixels
!        is less than MINSKY, an error code will result and
!        control will return to the main program.
!
! MAXSKY the maximum number of pixels allowed in the sky annulus.
!        This and the user's requested inner sky radius will later
!        determine the maximum permitted outer sky radius.
!
! MAXAP  the maximum number of star apertures allowed.
!
  integer, parameter :: precision = selected_real_kind(14)
  real(precision) :: maglim, magsq, wt, sumwt
  real(precision) :: apmag(maxap), area(maxap)
  real, dimension(:), allocatable :: sky
  real :: error(3), magerr(maxap)
  real :: par(maxap+2) = 0.0 
  real :: pmin(maxap+2) = tiny(0.0)
  real :: pmax(maxap+2) = huge(1.0)
  real :: lobad, skymod, skysig, skyskw, sigsq, skyvar, skymn, skymed
  real :: datum, r, rsq, fractn, edge, hibad, thresh, dum
  real :: phpadu, readns, xc, yc, dmag, watch, apmxsq
  real :: rinsq, rout, routsq, dysq
  integer :: i, j, k, l, naper, idum, nmag, lx, ly, nx, ny
  integer :: istar, mx, my, nsky, istat, nl, maxsky

  character(len=10) :: table = 'mphoto.opt' 
  character(len=26) :: lbl(maxap+2) = (/ &
              ' A1  RADIUS OF APERTURE  1',  & 
              ' A2  RADIUS OF APERTURE  2',   &
              ' A3  RADIUS OF APERTURE  3',   &
              ' A4  RADIUS OF APERTURE  4',   &
              ' A5  RADIUS OF APERTURE  5',   &
              ' A6  RADIUS OF APERTURE  6',   &
              ' A7  RADIUS OF APERTURE  7',   &
              ' A8  RADIUS OF APERTURE  8',   &
              ' A9  RADIUS OF APERTURE  9',   &
              ' AA  RADIUS OF APERTURE 10',   &
              ' AB  RADIUS OF APERTURE 11',   &
              ' AC  RADIUS OF APERTURE 12',   &
              ' IS       INNER SKY RADIUS',   &
              ' OS       OUTER SKY RADIUS'/)
!
!-----------------------------------------------------------------------
!
! SECTION 1
!
! Ascertain the name of the aperture photometry parameter table, and
! read it in.  Then set up all necessary variables for the forthcoming
! reductions. Finally, identify and open the input and output files.
!
  l = maxap + 1
  pmin(l) = 1.0

  call getopt(table, maxap+2, par, pmin, pmax, istat)
  if( istat /= 0 )then
     write(*,*) 'Error during read apertures radii. Using following values: '
     do i = 1, maxap + 2
        write(*,"(a,f10.2)") lbl(i),par(i)
     enddo
  endif
!
! Count up the number of apertures that will be used.  The first zero or
! negative number encountered terminates the list.
!
  naper = maxap
  apmxsq = -1.0
  do  i = 1, maxap
     if (par(i) <= 0.0 ) then
        naper = i - 1
        exit
     endif
     apmxsq = max(apmxsq, (par(i)+0.5)**2)
  enddo
!
! sky buffer allocation
!
  maxsky = int(max(Pi*((par(maxap+2)+1)**2 - par(maxap+1)**2),1.5))
  allocate(sky(maxsky))

!
! NAPER   is the number of apertures, whose radii are stored in
!         elements 1 through NAPER of the array PAR.
!
! APMXSQ  is the outermost edge of the largest aperture-- if the
!         distance squared of the center of a pixel from the centroid of
!         the star is greater than APMXSQ, then we know that no part
!         of the pixel is to be included in any aperture.
!
! Now define the other variables whose values are in the table.
!
  rinsq = max(par(maxap + 1), 0.0)**2      ! Inner sky radius squared
  routsq = maxsky/Pi + rinsq
  dum = par(maxap + 2)**2
  if (dum > routsq) then
     call stupid('   *** You have specified too big a sky annulus. ***')
     write (*,"(F10.2,A)") sqrt(routsq), &
          ' pixels is the largest outer sky radius currently permitted.' 
     return
  else if (dum <= rinsq) then
     call stupid ('Your outer sky radius is no bigger than the inner radius.')
     write (*,*) ' Please try again.'
     return
  else
     rout = par(maxap + 2)
     routsq = dum
  end if
!
! Inquire the name of the input data file with the stellar positions,
! and open it.
!
  call infile (2, coofil, istat)
  if (istat /= 0) then
     call stupid ('Error opening input file '//coofil)
     return
  end if
  call rdhead(2, nl, idum, idum, lobad, hibad, thresh, dum, phpadu, readns, dum)
  if (nl < 1) nl = 1
!
! Inquire file name for output aperture photometry results, and open
! the new file.
!             
  call outfil (3, magfil, istat)
  if (istat /= 0) then
     call stupid ('Error opening output file '//magfil)
     return
  end if
  call wrhead (3, 2, ncol, nrow, 6, lobad, hibad, thresh, par(1), phpadu, readns, 0.)
  readns = readns**2
!
! If progress is being monitored, type out column headers.
!
  if (watch > 0.5) &
       write (*,"(/13X, 'STAR', 5X, 'X', 7X, 'Y', 9X, 'MAG.(1)', 8X, 'SKY')")
!
! Initialize variables for the computation of the magnitude limit.
!
  maglim = 0.0
  magsq = 0.0
  sumwt = 0.0
  nmag = 0
!
!-----------------------------------------------------------------------
!
! SECTION 2
!
! Derive aperture photometry object by object.
!
! Get the coordinates of next object to be measured.
!
  lx = 1
  ly = 1
  nx = ncol
  ny = nrow

  do
     call rdstar (2, nl, istar, xc, yc, dmag, dum)
     if (istar < 0) exit
!      if (istar == 0) go to 2000
!
! Compute the limits of the submatrix.
!
     lx = max(1, int(xc-rout)+1)
     mx = min(ncol, int(xc+rout))
     ly = max(1, int(yc-rout)+1)
     my = min(nrow, int(yc+rout))
     edge = min(xc-0.5, (ncol+0.5)-xc, yc-0.5, (nrow+0.5)-yc)
!
! EDGE is the distance of the star's centroid from the outermost
! extremum of the array.
!
! Initialize star counts and aperture area.
!
     do  i = 1, naper
        apmag(i) = 0.0
!
! If this star aperture extends outside the array, the magnitude
! in this aperture will be no good.
!
        if (edge < par(i)) apmag(i) = - huge(1.0)            ! Null magnitude
        area(i) = 0.0
     enddo
!
! Now read through the submatrix, picking out the data we want.
!
     nsky = 0
!
     do  j = ly, my
        dysq = (j - yc)**2
!
        do i = lx,mx
           rsq = dysq + (i - xc)**2
           datum = d(i,j)
!
! Is this pixel within the sky annulus?
!
           if ( .not.((rsq < rinsq) .or. (rsq > routsq) .or. &
                (nsky > maxsky) .or. (datum < lobad) .or. &
                (datum > hibad))) then 
              nsky = nsky + 1
              sky(nsky) = datum
           endif
!
! The inclusion of partial pixels inside the aperture is done as
! follows:  if the distance of the center of the current pixel from the
! centroid of the star [radius vector r(i,j)] is exactly equal to the
! radius of the aperture [R(k)], then one-half of the counts in the
! pixel are included.  If r(i,j) < R(k)-0.5, then the entire pixel is
! included, while if r(i,j) > R(k)+0.5, the pixel is wholly excluded.
! In between, viz. for  R(k)-0.5 < r(i,j) < R(k)+0.5, the fraction of
! the counts included varies linearly.  Therefore a circular aperture
! is approximated by an irregular (not even convex) polygon.
!
! If this pixel falls completely outside the LARGEST aperture, go on
! to the next pixel.  Notice that APMXSQ has actually been defined
! as (R(k)+0.5)**2 for the largest value of R(k), in accordance with
! the formula used for the partial pixels.
!
           if (rsq <= apmxsq) then
              r = sqrt(rsq) - 0.5
!
              do k = 1, naper
!
! if this pixel falls completely outside THIS aperture, go on to the
! next aperture.
!
                 if (r <= par(k)) then
                    fractn = max(0.0, min(1.0,par(k) - r))
!
! fractn is the fraction of the pixel that falls inside the
! (irregular) aperture.
!
! If the pixel is bad, set the total counts in this aperture to a number
! so negative that it will never be positive again.
!                                                          ! Null magnitude
                    if (datum < lobad .or. datum > hibad ) apmag(k) = - huge(1.0) 
                    apmag(k) = apmag(k) + fractn*datum
                    area(k) = area(k) + fractn
                 endif
              enddo 
           endif
        enddo ! i
!
   enddo ! j
!
! We have accumulated the brightnesses of individual sky pixels in the
! one-dimensional array SKY.  Pixels falling above or below the BAD
! limits have already been eliminated.  Now sort SKY to place the
! pixels in order of increasing brightness.
!
   if (nsky < minsky)  then
      call stupid ('There aren''t enough pixels in the sky annulus.')
      write (*,*) ' Are you sure your bad pixel thresholds are all right?'
      write (*,*) ' If so, then you need a larger outer sky radius.'
      call tblank
      call clfile (2)
      call clfile (3)
      return
   end if
!      call quick (sky, nsky, index)
!
! Obtain the mode, standard deviation, and skewness of the peak in the
! sky histogram.
!
!      call mmm (sky, nsky, hibad, dum, datum, skymod, skysig, skyskw)
   call robustmean1(sky, nsky, skymn, skymed, skymod, skysig, skyskw)
   
   skyvar=skysig**2
   sigsq=skyvar/nsky

!
! SKYMOD has units of (ADU/pixel), and SKYSIG is the pixel-to-pixel
! scatter of SKYMOD, in units of (ADU/pixel).  SKYVAR is the
! variance (square of the standard deviation) of the sky brightness,
! (ADU/pixel)**2, and SIGSQ is the square of the standard error of the
! mean sky brightness.
!
! Subtract the sky from the integrated brightnesses in the apertures,
! convert the results to magnitudes, and compute standard errors.
!
   do i = 1, naper
!
! If the modal sky value could not be determined, set the magnitude
! to 99.999.
!
      if (skysig < -0.5) go to 2210
      apmag(i) = apmag(i) - skymod*area(i)
!
! If the star + sky is fainter than the sky, or if the star aperture
! extends beyond the limits of the picture, or if there is a bad pixel
! in the star aperture, set the magnitude to 99.999.
!
      if (apmag(i) <= 0.0) go to 2210
      error(1) = area(i)*skyvar
      error(2) = apmag(i)/phpadu
      error(3) = sigsq*area(i)**2
!
! These variables ERRORn are the respective variances (squares of the
! mean errors) for: (1) random noise inside the star aperture, including
! readout noise and the degree of contamination by other stars in the
! neighborhood, as estimated by the scatter in the sky values (this
! standard error increases as the square root of the area of the
! aperture); (2) the Poisson statistics of the observed star brightness;
! (3) the uncertainty of the mean sky brightness (this standard error
! increases directly with the area of the aperture).
!
      magerr(i) = min(9.999d0, 1.0857*sqrt(error(1) + error(2) + error(3))/apmag(i))
      apmag(i) = 25.0 - 2.5*log10(apmag(i))
      if (apmag(i) > 99.999) go to 2210
      go to 2220
2210  continue
      apmag(i) = 99.999
      magerr(i) = 9.999
2220  continue
      enddo 
!
!  NOTICE: a nice example of goto's
!
! Write out the answers.
!
      if (watch > 0.5) then
         write (*,"(12X, I5, 2F8.2, F9.3, ' +-', F6.3, F8.1)") &
              istar, xc, yc, apmag(1), magerr(1), skymod
      endif
      write (3,"(/1X, I5, 14F9.3)") istar, xc, yc, (apmag(i), i=1,naper)
      write (3,"(4X, F9.3, 2F6.2, F8.3, 11F9.3)") skymod, min(99.99,skysig), &
             min(999.99, max(-99.99,skyskw)), (magerr(i), i=1,naper)

      if (apmag(1) < 99.0) then
         wt = (2.0/(2.0 - dmag))*(100.0/magerr(1))**2
         maglim = maglim + wt*(apmag(1) - dmag)
         magsq = magsq + wt*(apmag(1) - dmag)**2
         sumwt = sumwt + wt
         nmag = nmag + 1
      endif
   enddo
!
!-----------------------------------------------------------------------
!
! Normal return.
!
! Estimate magnitude limit, close files, and return.
!

deallocate(sky)

call clfile (3)
call clfile (2)
if (sumwt <= 0.0) return
maglim = maglim/sumwt
magsq = magsq/sumwt - maglim**2
magsq = sqrt(max(0.0d0, magsq))

write (*,"('    Estimated magnitude limit (Aperture 1): ', F4.1,' +-', F4.1, ' per star.')") maglim, magsq
!      call stupid (line)
return
!
!-----------------------------------------------------------------------
!
end
